############################################
# NOTE: The replication is conducted in the order that tables and figures are placed in paper.
############################################

############################################ ################################################
######################################## PACKAGES ##########################################
############################################ ################################################

pacman::p_load(dplyr, tidyr, arsenal, devtools, ggstatsplot, estimatr, knitr, kableExtra, 
               ggplot2, texreg, rio, modelsummary, readstata13, stargazer, ggrepel, broom,
               rdrobust, tidyverse, tm, wordcloud, stringi, hunspell, syuzhet, tidyr, gridExtra,
               SnowballC, reshape2, stringi, sjPlot, lmtest, sandwich, marginaleffects,
               ggrepel, haven, lubridate)
if (!dir.exists("plots")) {
  dir.create("plots")
}

##### FIGURE 1 #####
setwd("~/Dropbox/EJPR_replication")
merkel <- read.csv("data/Google Search.csv")
merkel$Date <- as.Date(merkel$Day, format = "%Y-%m-%d")
merkel$Search_Merkel <- merkel$Angela.Merkel...Germany.
speech_text <- "Merkel's speech."
quarantine_text <- "Merkel's quarantine."
nearest_date_speech <- merkel$Date[which.min(abs(as.Date("2020-03-18") - merkel$Date))]
nearest_date_quarantine <- merkel$Date[which.min(abs(as.Date("2020-03-22") - merkel$Date))]

fig1 = ggplot(merkel, aes(x = Date, y = Search_Merkel)) +
  geom_line() +
  geom_vline(xintercept = as.Date("2020-03-18"), color = "red", linetype = "dashed", size = 1) +
  geom_vline(xintercept = as.Date("2020-03-22"), color = "blue", linetype = "dashed", size = 1) +  # New line for Merkel quarantine
  theme_minimal() +
  labs(x = "Date", y = "Google Search for Merkel") +
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  geom_label_repel(
    data = subset(merkel, Date == nearest_date_speech),
    aes(x = Date, y = Search_Merkel, label = speech_text),
    size = 3, nudge_x = -3, nudge_y = 10,
    box.padding = 0.5, segment.color = "transparent",
    arrow = arrow(length = unit(0.02, "npc"), type = "closed")
  ) +
  geom_label_repel(
    data = subset(merkel, Date == nearest_date_quarantine),
    aes(x = Date, y = Search_Merkel, label = quarantine_text),
    size = 3, nudge_x = 3, nudge_y = 10,  # Adjust positioning as needed
    box.padding = 0.5, segment.color = "transparent",
    arrow = arrow(length = unit(0.02, "npc"), type = "closed")
  ) +
  theme_minimal()
fig1
ggsave(filename = "plots/fig1.png", plot = fig1, width = 8, height = 6, dpi = 300)

##### FIGURE 2 #####
df = read_rds("data/Merkel_dataset.rds")

png("plots/fig2.png", width = 6, height = 4, units = "in", res = 300)  # Adjust width and height in inches
rdplot(
  y = df$Trust_in_Gov,
  x = df$day,
  c = 17,
  p = 1,
  ci = 95,
  binselect = "qsmv",
  kernel = "uniform",
  weights = df$webal,
  title = "Merkel's Speech and Political Trust",
  y.label = "Average Trust in Government (0-1 Scale)",
  x.label = "Day"
) + theme_minimal()
dev.off()

##### FIGURE 3 #####

covariate_names <- list(
  "(Intercept)" = "(Intercept)",
  "Speech" = "Speech",
  "Speech2" = "Speech",
  "Speech3" = "Speech",
  "Alternative_Speech1" = "Speech",
  "Alternative_Speech2" = "Speech",
  "Alternative_Speech3" = "Speech",
  'Party_AffiliationOther parties',
  "Speech:Party_AffiliationOther parties" = 'Speech x Other Pol.Parties',
  "Party_AffiliationSPD" = 'SPD',
  "Party_AffiliationFDP" = 'FDP',
  "Party_Affiliationdie Grunen" = 'die Grunen',
  "Party_Affiliationdie Linke" = 'die Linke',
  "Party_AffiliationOther" = 'Other',
  "women" = "Women",
  "age" = "Age",
  "has_partner" = "Has Partner",
  "intermediary_school" = "Intermed. School",
  "university_applied" = "Univ. Applied",
  "university_abitur" = "Univ. Abitur",
  "other_education" = "Other Education",
  "completed_vocational_training" = "Comp.Voc. Training",
  "completed_academic_degree" = "Comp. Acad. Degree",
  "other_vocational_degree" = "Other Voc. Degree",
  "full_time_employed" = "Full-Time",
  "part_time_employed" = "Part-Time",
  "education_related_status" = "Educational Emp.",
  "other_status" = "Other Emp.",
  "retired" = "Retired",
  "daily_internet" = "Daily Internet",
  "N_household_members" = "N HH Members",
  "day" = "Day",
  "I(day^2)" = "Day^2"
)

baseline<-lm_robust(Trust_in_Gov~Speech, fixed_effects = ~ factor(state_19), cluster=cluster, data=df, `se_type` = "stata")
extended<-lm_robust(Trust_in_Gov~Speech+day+I(day^2), cluster=cluster, fixed_effects = ~ factor(state_19),
                    data=df, `se_type` = "stata")
full<-lm_robust(Trust_in_Gov~Speech+
                  women+ age+
                  has_partner+intermediary_school+university_applied+
                  university_abitur+other_education+
                  completed_vocational_training+completed_academic_degree+
                  other_vocational_degree+full_time_employed+
                  part_time_employed+education_related_status+
                  other_status+retired+daily_internet+N_household_members+day+I(day^2), cluster=cluster, fixed_effects = ~ factor(state_19),
                data=df, `se_type` = "stata")
balanced<-lm_robust(Trust_in_Gov~Speech+
                      women+ age+
                      has_partner+intermediary_school+university_applied+
                      university_abitur+other_education+
                      completed_vocational_training+completed_academic_degree+
                      other_vocational_degree+full_time_employed+
                      part_time_employed+education_related_status+
                      other_status+retired+daily_internet+N_household_members+day+I(day^2),
                    fixed_effects = ~ factor(state_19), cluster=cluster, data=df, 
                    weights = webal, 
                    `se_type` = "stata")

coef_baseline <- broom::tidy(baseline, conf.int = TRUE) %>%
  filter(term == "Speech") %>%
  mutate(Model = "Baseline")

coef_extended <- broom::tidy(extended, conf.int = TRUE) %>%
  filter(term == "Speech") %>%
  mutate(Model = "Extended")

coef_full <- broom::tidy(full, conf.int = TRUE) %>%
  filter(term == "Speech") %>%
  mutate(Model = "Full")

coef_balanced <- broom::tidy(balanced, conf.int = TRUE) %>%
  filter(term == "Speech") %>%
  mutate(Model = "Balanced")

speech_coefs <- bind_rows(coef_baseline, coef_extended, coef_full, coef_balanced)
speech_coefs$Model <- factor(speech_coefs$Model, levels = c("Baseline", "Extended", "Full", "Balanced"))

fig3 = ggplot(speech_coefs, aes(x = Model, y = estimate, ymin = conf.low, ymax = conf.high, linetype = Model)) +
  geom_pointrange(aes(color = Model), size = 0.3) +  # Smaller points, color-coded by model
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +  # Zero line
  theme_minimal() + 
  scale_linetype_manual(values = c("solid", "solid", "solid", "solid")) +  # Custom line types
  scale_color_manual(values = c("Baseline" = "black", "Extended" = "black", "Full" = "black", "Balanced" = "red")) +  # Color Balanced in red
  labs(
    x = 'Model',
    y = 'Coefficient Estimate',
    title = 'The Effect of Merkel’s Speech on Political Trust'
  ) + 
  theme(
    legend.position = "none",  # Remove legend
    plot.title = element_text(size = 12, color = "black"),  # Title not in bold
    axis.title.x = element_text(color = "black"),
    axis.title.y = element_text(color = "black"),
    axis.text = element_text(color = "black")
  ) +
  coord_flip()
fig3
ggsave(filename = "plots/fig3.png", plot = fig3, width = 8, height = 6, dpi = 300)

##### FIGURE 4 #####

# Please refer to EJPR_forest_replication.R for Figure 4.

##### FIGURE 5 #####

fig5 = ggplot(df, aes(x=dDatum)) + 
  geom_bar() +
  geom_vline(xintercept = as.numeric(as.Date("2020-03-18")),
             color =  "red",
             linetype = "dashed") +
  labs(title = "Interviews per day before and after treatment",
       subtitle = "Red line indicates Merkel's speech on March 18, 2020",
       y = "Participants",
       x = "Day")+  theme_minimal()
fig5
ggsave(filename = "plots/fig5.png", plot = fig5, width = 8, height = 6, dpi = 300)

##### TABLE 1 #####
df$gender <- as.factor(df$gender_19)
df$educ_school <- as.factor(df$educ_school_19)
df$educ_job <- as.factor(df$educ_job_19)
df$occupation <- as.factor(df$occupation_19)
df$marital_status <- as.factor(df$marital_status_19)

df_desc = df %>% dplyr::select(gender, age, educ_school, educ_job,marital_status,number_hh_members_19,occupation,internet_usage_19, Party_Affiliation, Speech)
stargazer(df_desc)
summary(tableby(Speech ~ ., data = df_desc), title = "Balance Test of Covariates by Treatment", text="latex")

##### TABLE 2 #####

# NOTE: The code above produces both tables but they are divided into two in Overleaf due to space constraints.

##### TABLE 3 #####

df_desc = df %>% dplyr::select(AK46034, AA46041a10, AA46021, AE46108, Speech)
summary(df_desc, title = "Summary Statistics of Key Variables by Speech Treatment", text=TRUE)
summary(df_desc, title = "Summary Statistics of Key Variables by Speech Treatment", text="latex")

##### TABLE 4 #####
texreg(list(baseline, extended, full, balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Merkel's Speech and Trust in Federal Government", stars = c(0.001, 0.01, 0.05, 0.1),
       custom.model.names = c("Baseline", "Extended", "Full", 'Balanced'),
       fontsize="scriptsize",
       custom.coef.map = covariate_names,
       label = "tab:estimates-trust", 
       custom.gof.rows = list("State FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Clustered SE" = c("Yes", "Yes", "Yes","Yes"),
                              "Time" = c("No", "Yes", "Yes", "Yes"),
                              "Weight" = c("No", "No", "No","Yes")))
                                                             

#### TABLE 5 ####

dfmini <- subset(df, dDatum< as.Date("2020-03-22"))

model1<-lm_robust(Trust_in_Gov~Speech2+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day, fixed_effects = ~ factor(state_19), data=dfmini, cluster=cluster,weights = webal, `se_type` = "stata")

texreg(list(model1), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Restricting Sample to Pre-Quarantine Period", stars = c(0.001, 0.01, 0.05, 0.1), 
       custom.model.names = c("Balanced"),
       #omit.coef = "age|educ|day|gender|marital|number|occupation|internet",
       fontsize="scriptsize",
       custom.coef.map = covariate_names,
       #omit.coef = "age|educ|day|gender|marital|number|occupation|internet|occupation",
       label = "tab:quarantine", custom.gof.rows = list("State FE" = c("Yes"), "Covariates" = c("Yes"), "Clustered SE" = c("Yes"),  "Weight" = c("Yes"), "Time" = c("Yes")))

### TABLE 6 ####

cdu<-lm_robust(CDU_Placement~Speech+
                 women+ age+
                 has_partner+intermediary_school+university_applied+
                 university_abitur+other_education+
                 completed_vocational_training+completed_academic_degree+
                 other_vocational_degree+full_time_employed+
                 part_time_employed+education_related_status+
                 other_status+retired+daily_internet+N_household_members+day+I(day^2),
               fixed_effects = ~ factor(state_19), cluster=cluster, data=df, 
               weights = webal, 
               `se_type` = "stata")
spd<-lm_robust(SPD_Placement~Speech+
                 women+ age+
                 has_partner+intermediary_school+university_applied+
                 university_abitur+other_education+
                 completed_vocational_training+completed_academic_degree+
                 other_vocational_degree+full_time_employed+
                 part_time_employed+education_related_status+
                 other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
fdp<-lm_robust(FDP_Placement~Speech+
                 women+ age+
                 has_partner+intermediary_school+university_applied+
                 university_abitur+other_education+
                 completed_vocational_training+completed_academic_degree+
                 other_vocational_degree+full_time_employed+
                 part_time_employed+education_related_status+
                 other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
grunen<-lm_robust(Grunen_Placement~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
linke<-lm_robust(Linke_Placement~Speech+
                   women+ age+
                   has_partner+intermediary_school+university_applied+
                   university_abitur+other_education+
                   completed_vocational_training+completed_academic_degree+
                   other_vocational_degree+full_time_employed+
                   part_time_employed+education_related_status+
                   other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
afd<-lm_robust(AfD_Placement~Speech+
                 women+ age+
                 has_partner+intermediary_school+university_applied+
                 university_abitur+other_education+
                 completed_vocational_training+completed_academic_degree+
                 other_vocational_degree+full_time_employed+
                 part_time_employed+education_related_status+
                 other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
self<-lm_robust(Self_Placement~Speech+
                  women+ age+
                  has_partner+intermediary_school+university_applied+
                  university_abitur+other_education+
                  completed_vocational_training+completed_academic_degree+
                  other_vocational_degree+full_time_employed+
                  part_time_employed+education_related_status+
                  other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(cdu, spd, fdp, grunen, linke, afd, self), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Merkel's Speech and Ideological Placement", custom.model.names = c('CDU', "SPD", "FDP", "Grunen", "Linke", "AfD", "Self"),
       #omit.coef = "age|educ|day|gender|marital|number|occupation|internet",
       fontsize="scriptsize",
       custom.coef.map = covariate_names, label = "tab:otherparties",
       custom.gof.rows = list("State FE" = c("Yes", 'Yes', "Yes", 'Yes', "Yes", 'Yes', 'Yes'),
                              "Covariates" = c("Yes", 'Yes', "Yes", 'Yes', "Yes", 'Yes', 'Yes'),
                              "Clustered SE" = c("Yes", 'Yes', "Yes", 'Yes', "Yes", 'Yes', 'Yes'),
                              "Weight" = c("Yes", 'Yes', "Yes", 'Yes', "Yes", 'Yes', 'Yes'),
                              "Time" = c("Yes", 'Yes', "Yes", 'Yes', "Yes", 'Yes', 'Yes')))
                  
### TABLE 7 ####

mod1<-lm_robust(Trust_in_Gov~Speech+Party_Affiliation+
                  women+ age+
                  has_partner+intermediary_school+university_applied+
                  university_abitur+other_education+
                  completed_vocational_training+completed_academic_degree+
                  other_vocational_degree+full_time_employed+
                  part_time_employed+education_related_status+
                  other_status+retired+daily_internet+N_household_members+day, fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(mod1), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Merkel’s Speech, Political Affiliation and Political Trust", 
       stars = c(0.001, 0.01, 0.05, 0.1), custom.model.names = c("Control"),
       fontsize="scriptsize",
       custom.coef.map = covariate_names, label = "tab:interaction-party",
       custom.gof.rows = list("State FE" = c("Yes"),
                              "Covariates" = c("Yes"),
                              "Clustered SE" = c("Yes"),
                              "Weight" = c("Yes"),
                              "Time" = c("Yes")))

#### FIGURE 6 ####
df <- df %>%
  mutate(Party_Affiliation = case_when(
    AA43043 == 1 ~ "CDU/CSU",   
    TRUE ~ "Other parties"))

df$Party_Affiliation <- factor(df$Party_Affiliation, levels = c("CDU/CSU", "Other parties"))

df$Speech_Assignment <- factor(df$Speech, levels = c(0, 1), labels = c("Control", "Treatment"))
mod2 <- lm(Trust_in_Gov ~ Speech_Assignment * Party_Affiliation + 
               women + age + has_partner + intermediary_school + university_applied + 
               university_abitur + other_education + completed_vocational_training + 
               completed_academic_degree + other_vocational_degree + full_time_employed + 
               part_time_employed + education_related_status + other_status + retired + 
               daily_internet + N_household_members + day, data = df)
marginal_effects_plot <- plot_model(mod2, type = "int", terms = c("Speech_Assignment", "Party_Affiliation"),
                                    title = "Marginal Effects of Merkel’s Speech by Party Affiliation") +
  theme_minimal() + 
  ylab("Trust in Government") +
  scale_color_manual(values = c("black", "red"), 
                     labels = c("CDU/CSU", "Other parties")) +
  scale_fill_manual(values = c("black", "red"), 
                    labels = c("CDU/CSU", "Other parties"))
marginal_effects_plot
ggsave("plots/fig6.png", plot = marginal_effects_plot, device = "png")

#### TABLE 8 #####

df$Alternative_Speech1[df$dDatum<"2020-03-07"] = 0
df$Alternative_Speech1[df$dDatum=="2020-03-07"] = 0
df$Alternative_Speech1[df$dDatum>"2020-03-07"] = 1

df$Alternative_Speech2[df$dDatum<"2020-03-10"] = 0
df$Alternative_Speech2[df$dDatum=="2020-03-10"] = 0
df$Alternative_Speech2[df$dDatum>"2020-03-10"] = 1

df$Alternative_Speech3[df$dDatum<"2020-03-20"] = 0
df$Alternative_Speech3[df$dDatum=="2020-03-20"] = 0
df$Alternative_Speech3[df$dDatum>"2020-03-20"] = 1

speechdaytreatment<-lm_robust(Trust_in_Gov~Speech2+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
speechdayexclude<-lm_robust(Trust_in_Gov~Speech3+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
lockdown<-lm_robust(Trust_in_Gov~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2)+Lockdown, fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

covid<-lm_robust(Trust_in_Gov~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2)+Cases, fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
altspeech1<-lm_robust(Trust_in_Gov~Alternative_Speech1+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
altspeech2<-lm_robust(Trust_in_Gov~Alternative_Speech2+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")
altspeech3<-lm_robust(Trust_in_Gov~Alternative_Speech3+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(speechdaytreatment, speechdayexclude, lockdown, covid, altspeech1, altspeech2, altspeech3), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Robustness Checks with the Day of Speech", stars = c(0.001, 0.01, 0.05, 0.1),
       custom.model.names = c("Treat. on Speech Day", "Exclude Speech Day", "Lockdown", 'Covid',
                              'F.Date 1', 'F.Date 2', 'F.Date 3'),
       fontsize="scriptsize",
       custom.coef.map = covariate_names,
       label = "tab:estimates-rc-speech",
       custom.gof.rows = list("State FE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Clustered SE" = c("Yes", "Yes", "Yes","Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Weight" = c("Yes", "Yes", "Yes","Yes", "Yes", "Yes", "Yes")))

#### TABLE 9 ####
df$Trust_in_Court <- (df$AK46031-1) / 6

model1<-lm_robust(Trust_in_Court~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model2<-lm_robust(Satisfaction_CDU~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model3<-lm_robust(Populism1~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model4<-lm_robust(Populism2~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

model5<-lm_robust(Populism3~Speech+
                    women+ age+
                    has_partner+intermediary_school+university_applied+
                    university_abitur+other_education+
                    completed_vocational_training+completed_academic_degree+
                    other_vocational_degree+full_time_employed+
                    part_time_employed+education_related_status+
                    other_status+retired+daily_internet+N_household_members+day+I(day^2), fixed_effects = ~ factor(state_19), cluster=cluster,weights = webal, data=df, `se_type` = "stata")

texreg(list(model1, model2, model3, model4, model5), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Speech Treatment and Other Political Variables",
       custom.model.names = c("TrustCourt", "CDU", "Pop1", "Pop2", 'Pop3'),
       omit.coef = "women|age|partner|school|applied|education|training|degree|time|emp|retired|internet|member|day|",
       fontsize="scriptsize",
       custom.coef.map = covariate_names, label = "tab:otherpolitical",
       custom.gof.rows = list("Covariates" = c("Yes", "Yes", "Yes", "Yes", "Yes"),
                              "State FE" = c("Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Clustered SE" = c("Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Weights" = c("Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("Yes", "Yes", "Yes", "Yes", "Yes"),
                              "Time^2" = c("Yes", "Yes", "Yes", "Yes", "Yes")))

#### FIGURE 7 #####

mobility_ger = import("data/changes-visitors-covid.csv")
mobility_ger = mobility_ger %>%
  filter(Entity =="Germany")
mobility_ger$Day = as.Date(mobility_ger$Day, "%Y-%m-%d")
mobility_ger = mobility_ger %>%
  filter(Day >= as.Date("2020-03-01")  & Day < as.Date("2020-04-01"))
mobility_ger = mobility_ger %>% rename("Residential"="residential",
                               "Retail" = "retail_and_recreation",
                               "Grocery" = "grocery_and_pharmacy",
                               "Transit" = "transit_stations",
                                "Parks" = "parks",
                               "Workplaces" = "workplaces")

df <- mobility_ger %>%
  select(Day, Residential, Retail, Grocery, Transit, Parks, Workplaces) %>%
  gather(key = "variable", value = "value", -Day)
head(df)

fig7= ggplot(df, aes(x = Day, y = value)) + 
  geom_line(aes(color = variable, linetype = variable)) + 
  geom_vline(xintercept = as.numeric(as.Date("2020-03-18")), color="red", linetype = "dashed") +
  theme_minimal()+
  scale_x_date(date_breaks = "1 day", date_labels = "%b %d") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  labs(y = "Mobility", x= "Date")+
  scale_color_manual(values = c("gray0", "gray15", "gray30",
                                "gray45", "gray60", "gray75"))
fig7
ggsave(filename = "plots/fig7.png", plot = fig7, width = 8, height = 6, dpi = 300)


#### TABLE 10 #####

mobility <- read.csv("data/2020_DE_Region_Mobility_Report.csv") %>%
  select(-country_region_code, -country_region, -sub_region_2, -metro_area, -census_fips_code) %>%
  mutate(date = as.Date(date, "%Y-%m-%d")) %>%
  filter(date >= as.Date("2020-03-01") & date < as.Date("2020-04-01")) %>%
  rename(
    residential = residential_percent_change_from_baseline,
    retail_and_recreation = retail_and_recreation_percent_change_from_baseline,
    grocery_and_pharmacy = grocery_and_pharmacy_percent_change_from_baseline,
    transit_stations = transit_stations_percent_change_from_baseline,
    parks = parks_percent_change_from_baseline,
    workplaces = workplaces_percent_change_from_baseline
  )

df <- mobility %>%
  select(date, sub_region_1, retail_and_recreation, grocery_and_pharmacy, residential, transit_stations, parks, workplaces) %>%
  pivot_longer(cols = -c(date, sub_region_1), names_to = "variable", values_to = "value")

mobility <- mobility %>%
  mutate(
    Speech = ifelse(date > as.Date("2020-03-18"), 1, 0),
    cluster = as.factor(paste(date, sub_region_1, sep = "")),
    dayspassed = as.numeric(date - as.Date("2020-03-01")),
    Lockdown = ifelse(date >= as.Date("2020-03-22"), 1, 0)
  )

cases_by_date <- c(
  "2020-03-01" = 1619, "2020-03-02" = 2039, "2020-03-03" = 2555, "2020-03-04" = 3562, "2020-03-05" = 5446,
  "2020-03-06" = 7677, "2020-03-07" = 33848, "2020-03-08" = 11623, "2020-03-09" = 12894, "2020-03-10" = 17069,
  "2020-03-11" = 24205, "2020-03-12" = 33297, "2020-03-13" = 45172, "2020-03-14" = 62757, "2020-03-15" = 78386,
  "2020-03-16" = 90141, "2020-03-17" = 114586, "2020-03-18" = 151026, "2020-03-19" = 194267, "2020-03-20" = 242834,
  "2020-03-21" = 291484, "2020-03-22" = 331739, "2020-03-23" = 358907, "2020-03-24" = 403671, "2020-03-25" = 462182,
  "2020-03-26" = 530324, "2020-03-27" = 601093, "2020-03-28" = 672833, "2020-03-29" = 729448, "2020-03-30" = 766392,
  "2020-03-31" = 815931
)
mobility$Cases <- as.numeric(cases_by_date[as.character(mobility$date)])

# Assign East/West regions based on sub_region_1
west_states <- c("Baden-Württemberg", "Bavaria", "Bremen", "Hamburg", "Hessen",
                 "Lower Saxony", "North Rhine-Westphalia", "Rhineland-Palatinate",
                 "Saarland", "Schleswig-Holstein")
east_states <- c("Berlin", "Brandenburg", "Mecklenburg-Vorpommern", "Saxony",
                 "Saxony-Anhalt", "Thuringia")

mobility <- mobility %>%
  mutate(
    EastWest = case_when(
      sub_region_1 %in% west_states ~ "West",
      sub_region_1 %in% east_states ~ "East",
      TRUE ~ NA_character_
    )
  )

range01 <- function(x){(x-min(x))/(max(x)-min(x))}

mobility_st <- 
  mobility %>% 
  mutate(grocery_and_pharmacy_s = range01(grocery_and_pharmacy),
         parks_s = range01(parks),                               
         residential_s = range01(residential),
         retail_and_recreation_s = range01(retail_and_recreation),
         transit_stations_s = range01(transit_stations),
         workplaces_s = range01(workplaces))

# Compliance at State level
model1<-lm_robust(grocery_and_pharmacy_s~Speech+dayspassed+I(dayspassed^2)+Lockdown+Cases, fixed_effects = ~ factor(sub_region_1), data=mobility_st,cluster=date, `se_type` = "stata")
model2<-lm_robust(parks_s~Speech+dayspassed+I(dayspassed^2)+Lockdown+Cases, fixed_effects = ~ factor(sub_region_1), data=mobility_st,cluster=date, `se_type` = "stata")
model3<-lm_robust(residential_s~Speech+dayspassed+I(dayspassed^2)+Lockdown+Cases, fixed_effects = ~ factor(sub_region_1), data=mobility_st,cluster=date, `se_type` = "stata")
model4<-lm_robust(retail_and_recreation_s~Speech+dayspassed+I(dayspassed^2)+Lockdown+Cases, fixed_effects = ~ factor(sub_region_1), data=mobility_st,cluster=date, `se_type` = "stata")
model5<-lm_robust(transit_stations_s~Speech+dayspassed+I(dayspassed^2)+Lockdown+Cases, fixed_effects = ~ factor(sub_region_1), data=mobility_st,cluster=date, `se_type` = "stata")
model6<-lm_robust(workplaces_s~Speech+dayspassed+I(dayspassed^2)+Lockdown+Cases, fixed_effects = ~ factor(sub_region_1), data=mobility_st,cluster=date, `se_type` = "stata")

texreg(list(model1, model2, model3, model4, model5, model6), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Merkel's Speech and Mobility at State-Date Level", stars = c(0.001, 0.01, 0.05, 0.1),
       label = "tab:estimates-mobility", custom.gof.rows = list("State FE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), "Clustered SE" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), "Time" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), "Time^2" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))


#### FIGURE 8 ####
rutte <- read.csv("data/multiTimeline_nl.csv", sep = ",", header = TRUE)
rutte$Date <- as.Date(rutte$Day, format = "%Y-%m-%d")
rutte$Search_rutte <- rutte$Rutte
speech_text <- "Rutte's speech."
nearest_date_speech <- rutte$Date[which.min(abs(as.Date("2020-03-16") - rutte$Date))]

fig8 = ggplot(rutte, aes(x = Date, y = Search_rutte)) +
  geom_line() +
  geom_vline(xintercept = as.Date("2020-03-16"), color = "red", linetype = "dashed", size = 1) +
  theme_minimal() +
  labs(x = "Date", y = "Google Search for Rutte") +
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  geom_label_repel(
    data = subset(rutte, Date == nearest_date_speech),
    aes(x = Date, y = Search_rutte, label = speech_text),
    size = 3, nudge_x = -3, nudge_y = 10,
    box.padding = 0.5, segment.color = "transparent",
    arrow = arrow(length = unit(0.02, "npc"), type = "closed")
  ) +
  theme_minimal()
fig8
ggsave(filename = "plots/fig8.png", plot = fig8, width = 8, height = 6, dpi = 300)

#### FIGURE 9 ####
politics <- read_dta("data/cv20l_EN_1.0p.dta")
back <- read_dta("data/avars_202003_EN_1.0p.dta") 
dat <- left_join(politics, back, by = "nomem_encr")
rm(politics, back)

dat$gov_trust <- dat$cv20l013 %>%
  na_if(-9) %>%
  na_if(-14) %>%
  na_if(-15)

dat$gov_trust <- (dat$gov_trust - min(dat$gov_trust, na.rm = TRUE)) / 
  (max(dat$gov_trust, na.rm = TRUE) - min(dat$gov_trust, na.rm = TRUE))

dat <- dat %>%
  filter(
    cv20l_m3 == 202003,   
    !is.na(cv20l_m3),
    cv20l298 != " "     
  ) %>%
  
  mutate(
    date = dmy(cv20l298),               
    speech = ifelse(date > as.Date("2020-03-16"), 1, 0),  
    day = as.numeric(date - min(date, na.rm = TRUE))      
  ) %>%
  select(nomem_encr,gov_trust, date, speech,day)


png("plots/fig9.png", width = 6, height = 4, units = "in", res = 300)  # Adjust width and height in inches
rdplot(y = dat$gov_trust,
       x = dat$day,
       c = 14,
       p = 1,
       ci= 95,
       binselect = "qsmv",
       kernel = "uniform",
       title = "Rutte's Speech and Political Trust",
       y.label = "Average Trust in Government (0-1 Scale)",
       x.label = "Day")

#### FIGURE 10 ####
johnson <- read.csv("data/multiTimeline_uk.csv", sep = ",", header = TRUE)
johnson$Date <- as.Date(johnson$Day, format = "%Y-%m-%d")
johnson$Search_johnson <- johnson$Boris.Johnson...United.Kingdom.
speech_text <- "Johnson's speech."
quarantine_text <- "Johnson tests positive."
nearest_date_speech <- johnson$Date[which.min(abs(as.Date("2020-03-23") - johnson$Date))]
nearest_date_quarantine <- johnson$Date[which.min(abs(as.Date("2020-03-27") - johnson$Date))]

fig10 <- ggplot(johnson, aes(x = Date, y = Search_johnson)) +
  geom_line() +
  geom_vline(xintercept = as.Date("2020-03-23"), color = "red", linetype = "dashed", size = 1) +
  geom_vline(xintercept = as.Date("2020-03-27"), color = "blue", linetype = "dashed", size = 1) +
  theme_minimal() +
  labs(x = "Date", y = "Google Search for Johnson") +
  scale_x_date(date_labels = "%b %d", date_breaks = "1 week") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  geom_label_repel(
    data = subset(johnson, Date == nearest_date_speech),
    aes(x = Date, y = Search_johnson, label = speech_text),
    size = 3, nudge_x = -3, nudge_y = 10,
    box.padding = 0.5, segment.color = "transparent",
    arrow = arrow(length = unit(0.02, "npc"), type = "closed")
  ) +
  geom_label_repel(
    data = subset(johnson, Date == nearest_date_quarantine),
    aes(x = Date, y = Search_johnson, label = quarantine_text),
    size = 3, nudge_x = 3, nudge_y = 10,
    box.padding = 0.5, segment.color = "transparent",
    arrow = arrow(length = unit(0.02, "npc"), type = "closed")
  )

fig10
ggsave(filename = "plots/fig10.png", plot = fig10, width = 8, height = 6, dpi = 300)

#### FIGURE 11 ####
df20192020 <- read_dta("data/apsp_a19m20_eul_pwta22.dta")
df20192020$Date <- as.Date(df20192020$REFDTE, format="%d%m%Y")
df_march <- subset(df20192020, format(Date, "%Y") == "2020" & format(Date, "%m") == "03")
df_march[df_march == -9 | df_march == -8] <- NA
df_march$days <- as.numeric(format(df_march$Date, "%d"))
df_march$Treatment_PM <- ifelse(df_march$Date > as.Date("2020-03-23"), 1, 0)

df_march$SATIS_norm <- (df_march$SATIS - min(df_march$SATIS, na.rm = TRUE)) / 
  (max(df_march$SATIS, na.rm = TRUE) - min(df_march$SATIS, na.rm = TRUE))

png("plots/fig11.png", width = 6, height = 4, units = "in", res = 300)  # Adjust width and height in inches
rdplot(
  y = df_march$SATIS_norm,
  x = df_march$days,
  c = 22,
  p = 1,
  ci = 95,
  binselect = "qsmv",
  kernel = "uniform",
  title = "RD Plot for Johnson's Speech",
  y.label = "Average Life Satisfaction (0-1)",
  x.label = "Day"
) + theme_minimal()
dev.off()

#### FIGURE 12 #####

speech_files <- list("speech1.txt", "speech2.txt", "speech3.txt")
speakers <- c("Rutte", "Merkel", "Johnson")
additional_stopwords <- c("can", "will", "must", "also", "even", "like", "one", 'many', 'much')

preprocess_text <- function(text, additional_stopwords) {
  corpus <- Corpus(VectorSource(text))
  corpus <- tm_map(corpus, content_transformer(tolower))
  corpus <- tm_map(corpus, removePunctuation)
  corpus <- tm_map(corpus, removeNumbers)
  corpus <- tm_map(corpus, removeWords, c(stopwords("english"), additional_stopwords))
  corpus <- tm_map(corpus, stripWhitespace)
  corpus <- tm_map(corpus, content_transformer(stri_trim_both))
  
  clean_text <- sapply(corpus, function(x) {
    words <- unlist(strsplit(as.character(x), " "))
    correct_words <- words[hunspell_check(words)]
    paste(correct_words, collapse = " ")
  })
  return(Corpus(VectorSource(clean_text)))
}

sentiments <- list()
for (i in seq_along(speech_files)) {
  speech_text <- paste(readLines(speech_files[[i]], warn = FALSE), collapse = " ")
  corpus <- preprocess_text(speech_text, additional_stopwords)
  sentiment_scores <- get_nrc_sentiment(as.character(corpus[[1]]))
  sentiments[[i]] <- colSums(sentiment_scores)
}

sentiments_df <- do.call(rbind, sentiments)
sentiments_df <- as.data.frame(sentiments_df)
sentiments_df$Speaker <- speakers

sentiments_long <- sentiments_df %>%
  pivot_longer(cols = -Speaker, names_to = "Sentiment", values_to = "Count") %>%
  group_by(Speaker) %>%
  mutate(Proportion = Count / sum(Count))

sentiments_long$Speaker <- factor(sentiments_long$Speaker, levels = c("Merkel", "Rutte", "Johnson"))

custom_colors <- c("Merkel" = "red", "Rutte" = "blue", "Johnson" = "black")

fig12_1 <- ggplot(sentiments_long, aes(x = Sentiment, y = Count, fill = Speaker)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = custom_colors) +
  theme_minimal() +
  labs(title = "Sentiment Count per Speaker", x = "Sentiment", y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

fig12_2 <- ggplot(sentiments_long, aes(x = Sentiment, y = Proportion, fill = Speaker)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = custom_colors) +
  theme_minimal() +
  labs(title = "Sentiment Proportion per Speaker", x = "Sentiment", y = "Proportion") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

fig12 <- grid.arrange(fig12_1, fig12_2, ncol = 1)

# Save the combined plot as a single image file
ggsave("plots/fig12.png", fig12, width = 10, height = 12, dpi = 300)

#### FIGURE 13 ####
speech_files <- c("speech1.txt", "speech2.txt", "speech3.txt")
speakers <- c("Rutte's Speech", "Merkel's Speech", "Johnson's Speech")

speeches <- data.frame(
  text = sapply(speech_files, function(file) paste(readLines(file), collapse = " ")),
  speaker = speakers,
  stringsAsFactors = FALSE
)

additional_stopwords <- c("can", "also", "even", "like", "one", "many", "much")

anxiety_terms <- c(
  "anxiety", "fear", "panic", "worry", "distress", "nervous", "tension", "uncertainty", "concern", 
  "apprehension", "dread", "alarm", "fright", "unease", "restless", "agitation", "anguish", 
  "hesitation", "insecurity", "turmoil", "crisis", "stress", "paranoia", "doubt", "wary", "fearful", 
  "trepidation", "foreboding", "disquiet", "discomfort", "skepticism", "hesitance", "pressure", 
  "nervousness", "worries", "apprehensive", "disturbed", "alarmist", "upset", "disconcert", 
  "tense", "fraught", "irritability", "jitters", "worrisome", "fret", "disquieting", "alarming", 
  "concerned", "disturbance", "unsettling", "stressed", "suspicion", "anguished", "mistrust",
  "threat", "danger", "afraid", "misfortune", "misery"
)

calmness_terms <- c(
  "calm", "peace", "assure", "reassure", "reassuring", "composed", "relaxed", "soothe", "steady", "tranquil", 
  "comfort", "serenity", "equanimity", "confidence", "confident", "security", "relief", "stable", "stillness", 
  "balanced", "poised", "ease", "harmony", "self-assured", "placid", "quiet", 
  "content", "gentle", "patience", "collected", "certainty", "tranquility", "relaxation", "control", 
  "assured", "untroubled", "settled", "calming", "repose", "cool-headed", "undisturbed", "comforted", 
  "secure", "grounded", "centered", "relieving", "mindful", "quieting", "placidity", "trust",
  "harmonious", "mild", "serene", "protect", "smooth", "halcyon", "hushed", "pacific", "bland", "cool", "safe"
)

accommodation_terms <- c(
  "citizens", "public", "community", "people", "residents", "society", "families", "individuals",
  "workers", "children", "elderly", "youth", "parents", "friends", "neighbors", "partners",
  "colleagues", "relatives", "households", "population", "inhabitants", "crowd", "audience",
  "everyone", "volunteers", "helpers", "supporters", "beneficiaries", "participants", "staff",
  "students", "patients", "recipients", "all", "neighbors", "teams", "collective", "cooperation",
  "inclusion", "us", "together", "each other", "mutual", "societal", "group", "public health",
  "healthcare workers", "carers", "frontline", "helpers", "assistance", "unity", "shared", "solidarity",
  "voters", "employees", "workforce", "jobseekers", "enterprises", "entrepreneurs", "producers", 
  "taxpayer", "civilians", "Brits", "refugees", "asylum_seekers"
)

accessibility_terms <- c("simple", "clear", "easy", "concise", "readable", "accessible", "direct", "understandable")
alignment_terms <- c("agree", "consensus", "support", "align", "cooperate", "coordinate", "harmonize", "unite", 'solidarity')

preprocess_text <- function(text, additional_stopwords) {
  corpus <- Corpus(VectorSource(text))
  corpus <- tm_map(corpus, content_transformer(tolower))
  corpus <- tm_map(corpus, removePunctuation)
  corpus <- tm_map(corpus, removeNumbers)
  corpus <- tm_map(corpus, removeWords, c(stopwords("english"), additional_stopwords))
  corpus <- tm_map(corpus, stripWhitespace)
  return(corpus)
}

create_tdm <- function(speech, additional_stopwords) {
  corpus <- preprocess_text(speech, additional_stopwords)
  tdm <- TermDocumentMatrix(corpus)
  as.matrix(tdm)
}

calculate_measure_score <- function(tdm, terms) {
  if (is.null(dim(tdm))) return(0)
  word_frequencies <- rowSums(tdm)
  
  score <- 0
  for (term in terms) {
    pattern <- paste0("\\b", term, ".*\\b")  # Match any word starting with the term
    matched_words <- grep(pattern, names(word_frequencies), value = TRUE)
    score <- score + sum(word_frequencies[matched_words])
  }
  return(score)
}

measure_scores <- data.frame(Speech = speeches$speaker)
for (i in 1:nrow(speeches)) {
  tdm <- create_tdm(speeches$text[i], additional_stopwords)
  word_count <- length(unlist(strsplit(speeches$text[i], "\\s+")))
  
  measure_scores[i, "Accessibility"] <- calculate_measure_score(tdm, accessibility_terms) / word_count
  measure_scores[i, "Allay"] <- (calculate_measure_score(tdm, calmness_terms) - 
                                   calculate_measure_score(tdm, anxiety_terms)) / word_count
  measure_scores[i, "Accommodation"] <- calculate_measure_score(tdm, accommodation_terms) / word_count
  measure_scores[i, "Alignment"] <- calculate_measure_score(tdm, alignment_terms) / word_count
}

measure_scores$Speech <- factor(measure_scores$Speech, levels = c("Merkel's Speech", "Rutte's Speech", "Johnson's Speech"))
measure_scores_long <- melt(measure_scores, id.vars = "Speech", variable.name = "Measure", value.name = "Score")

custom_colors <- c("red", "blue", "grey", "black")

fig13 <- ggplot(measure_scores_long, aes(x = Speech, y = Score, fill = Measure)) +
  geom_bar(stat = "identity", position = "dodge") +
  theme_minimal() +
  labs(
    title = "Scores by Measure in Speeches",
    x = "Speech",
    y = "Score",
    fill = "Category"
  ) +
  scale_fill_manual(values = custom_colors) +
  theme(legend.position = "right") +
  facet_wrap(~ Measure, scales = "free_y")
fig13
# Save the plot as fig13.png
ggsave("plots/fig13.png", fig13, width = 10, height = 6, dpi = 300)
